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ABSTRACT 

The molecular evolution that occurs in collapsing prestellar cores is investigated. To model 
the dynamics, we adopt the Larson-Penston solution and analogues with slower rates of collapse. 
For the chemistry, we utilize the new standard model (NSM) with the addition of deuterium 
fractionation and grain-surface reactions treated via the modified rate approach. The use of 
surface reactions distinguishes the present work from our previous model. We find that these 
reactions efficiently produce H2O, H2CO, CH3OH, N2, and NH3 ices. In addition, the surface 
chemistry influences the gas-phase abundances in a variety of ways. For example, formation of 
molecular nitrogen on grain surfaces followed by desorption into the gas enhances the abundance 
of this gas-phase species and its daughter products N 2 H + and NH 3 . The current reaction network 
along with the Larson-Penston solution allows us to reproduce satisfactorily most of the molecular 
column densities and their radial distributions observed in L1544. The agreement tends to 
worsen with models that include strongly delayed collapse rates. Inferred radial distributions in 
terms of fractional abundances are somewhat harder to reproduce. In addition to our standard 
chemical model, we have also run a model with the UMIST gas-phase chemical network. The 
abundances of gas-phase sulphur-bearing molecules such as CS and CCS are significantly affected 
by uncertainties in the gas-phase chemical network. In all of our models, the column density of 
N2H 4 " monotonically increases as the central density of the core increases during collapse from 
3x 10 4 cm~ 3 to 3x 10 7 cm -3 . Thus, the abundance of this ion can be a probe of evolutionary stage. 
Molecular D/H ratios in assorted cores are best reproduced in the Larson-Penston picture with 
the conventional rate coefficients for fractionation reactions. If we adopt the newly measured and 
calculated rate coefficients, the D/H ratios, especially N2D + /N2H + , become significantly lower 
than the observed values. 

Subject headings: stars: formation — ISM: molecules — ISM: clouds — ISM: individual (L1544) 

1. Introduction estimate the density and temperature distribution 

of such cores from dust continuum observations 
Studies of star formation are largely based on at m illimeter/submillimeter and infrared wave- 

the observation of molecular cloud cores. We can 
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lengths. Continuum observations with SCUBA 
and ISO have especially improved our knowledge 
of prestellar cores, which are cold dense cores with- 
out IRAS sources. Here the density distribution 
has been found to be rather flat in the inner re- 
gions and roughly consistent with a power law of 
the form r~ 2 at larger radii (Ward-Thompson et 
al. 1994; Andre, Ward-Thompson, & Barsony 
2000; Bacmann et al. 2000; Shirley et al. 2000; 
Evans et al. 2001). As well as being an indepen- 
dent probe of the density and temperature dis- 
tribution, molecular line emission can be a probe 
of the physical dynamics and chemical kinetics in 
prestellar cores. The interpretation of line obser- 
vations is more complicated than that of the con- 
tinuum, however, because molecular abundances 
may vary within a core and also with time. Recent 
observations have illustrated the different depen- 
dences of molecular emission on position in a vari- 
ety of prestellar cores. For example, the ion N 2 H+ 
is relatively compact and shows a well-defined 
peak while CO, CS, and CCS are more diffuse and 
sometimes show a central hole. Such emission pat- 
terns indicate that the latter species are depleted 
at the center while N2H" 1 " is not (Kuiper, Langer, 
& Velusamy 1996; Willacy, Langer, & Velusamy 
1998; Casein et al. 1999; Ohashi et al. 1999; Lee, 
Myers, & Tafalla 2001). Moreover, strong deu- 
terium fractionation, an indicator of the depletion 
of heavy molecules, has also been detected (Hirota 
et al. 1998; Caselli et al. 1999). An understand- 
ing of the chemistry occurring in prestellar cores 
would appear to be essential to develop a detailed 
knowledge of the conditions in the early stages of 
star formation. 

Theoretical studies of the chemistry in prestel- 
lar cores have been performed by several groups. 
Bergin & Langer (1997) calculated molecular evo- 
lution in a fluid element in which the density 
increases as a function of time. These authors 
included gas-phase reactions and adsorption and 
desorption of molecules on grain surfaces. In their 
model, species such as CO are depleted from the 
gas due to adsorption onto grains, but N 2 H + sur- 
vives because of the lower binding energy of its 
parent species N 2 on grain surfaces compared with 
other heavy molecules. Aikawa et al. (2001) (here- 
after Paper I) utilized a similar chemical model 
but with multiple fluid elements collapsing accord- 
ing to the semi-analytical Larson-Penston solution 



(Larson 1969; Penston 1969) in order to obtain 
the spatial distributions of molecules as a function 
of time. The column densities of CO, CS, and CCS 
were found to have a central hole, while N 2 H + 
was found to be centrally peaked. Paper I also 
investigated molecular column densities in more 
slowly collapsing cores than the Larson-Penston 
flow to mimic the effect of magnetic support or 
rotation. It was found that the CCS abundance 
is significantly lower in these cores, and that a 
Larson-Penston core best reproduces the observed 
molecular abundances in L1544. Since molecular 
abundances are determined by a balance between 
chemical and collapse time scales, the chemistry 
represents a probe of the collapse time scale and 
thus of the physical mechanism of star formation. 

Li ct al. (2002) went one step further, directly 
coupling a chemical model of prestellar cores with 
a magneto-hydrodynamic simulation of the dy- 
namics. The results of Li et al. (2002) show 
that a core model with magnetic support repro- 
duces the observed molecular column densities in 
L1544 better than a non-magnetized model. The 
reasons for the different conclusions — L1544 is 
best reproduced by the Larson-Penston model in 
our Paper I but by the magnetized model in Li et 
al. (2002) — are, however, not fully known. The 
analysis of Li et al. (2002) actually has two advan- 
tages: inclusion of grain-surface reactions and di- 
rect coupling between the chemistry and magneto- 
hydrodynamic simulation via ambipolar diffusion. 
But, in order to include these two processes, Li 
ct al. (2002) were forced to reduce the chemical 
network to save computational time. 

This paper is an update of the earlier anal- 
ysis of Paper I. Here we investigate molecular 
distributions in a core collapsing according to 
the Larson-Penston solution and its analogues 
with delayed collapse using the new standard 
model (NSM) (Tcrzieva & Herbst 1998) with deu- 
terium fractionation, and grain-surface reactions 
that are treated by the modified rate approach 
(Stantcheva, Caselli, & Herbst 2001). With the 
semi-analytical Larson-Penston solution as a phys- 
ical model, we can easily incorporate the extensive 
chemical network needed. There are three goals 
of this paper. Firstly, we hope to clarify what 
is the main cause of the different conclusions be- 
tween our previous work and Li et al. (2002). If 
the main cause is the inclusion of grain-surface 
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reactions, our present results will become simi- 
lar to those of Li et al. (2002). Secondly, we 
will describe the detailed distribution of gaseous 
and adsorbed molecules in assorted evolutionary 
stages of prestellar cores in order to determine 
their value as probes of core evolution. Deuter- 
ated species are of special interest because the 
D /H ratio has been used as a probe of the degree 
of ionization (Caselli et al. 1998) and the deple- 
tion of heavy gas- phase molecules (Brown & Millar 
1989; Roberts & Millar 2000). The distribution of 
adsorbed species in prestellar cores has not been 
observed, but it provides an important initial con- 
dition for warm protostellar cores, in which icy 
mantles will be desorbed to start a second cycle 
of gas-phase reactions. Thirdly, we will compare 
our results with observations of several prestel- 
lar cores including L1544 and less dense prestellar 
cores such as L1521E, which are often described as 
'starless cores' rather than 'prestellar cores'. Con- 
sidering the uncertainties both in observational 
data and chemical predictions, which will be illus- 
trated below, we need to compare observational 
data and theoretical models for as many sources 
as feasible. In addition, comparison of various 
prestellar cores with model results should serve to 
tell us about the formation and evolution of cores. 

The remainder of the paper is organized as fol- 
lows. The description of our dynamical and chem- 
ical models is given in §2. In §3, we describe 
the evolution and distribution of molecular abun- 
dances in collapsing cores. We compare our re- 
sults with those of our previous paper, in which 
grain-surface reactions were not included. The 
dependence of our results on collapse time scale, 
sticking probability, initial gas density, and gas- 
phase chemical network is also investigated. In 
§4, our model results are compared with observa- 
tional data, as well as with the results of Li et 
al. (2002), while a summary is contained in §5. 

2. Model 

The distributions of density and velocity of the 
core, as obtained from the Larson-Penston model, 
are shown in Figure 1 of Paper I at several evo- 
lutionary stages. Because the Larson-Penston so- 
lution is self-similar, the "shape" of the density 
distribution is always the same: it is almost flat 
at smaller radii, and follows a power-law n oc r~ 2 



at larger radii. The initial central density of the 
core, which is nearly homogeneous at radii through 
10000 AU, is assumed to be n H = 2 x 10 4 cm" 3 
(i.e., n(H2)= 1 x 10 4 cm -3 ). The core initially 
extends to the radius of 3.2 x 10 4 AU, at which 
the density is jih = 5.3 x 10 3 cm~ 3 . The subse- 
quent collapse of the core occurs in almost a free- 
fall time scale, with the central density increasing 
to riH = 3 x 10 7 cm~ 3 and the outermost radius de- 
creasing to 1.36 x 10 4 AU in 2.00 x 10 5 yr. During 
the evolution, the core is isothermal with T = 10 
K. The Larson-Penston flow is an asymptotic so- 
lution of a non-magnetized collapsing core; hydro- 
dynamic simulations and normal mode analysis 
show that the flow asymptotically converges to the 
Larson-Penston flow (Foster & Chevalier 1993; 
Hanawa & Nakayama 1997). Although it is a good 
approximation of collapsing dense cores especially 
in their inner regions, the Larson-Penston flow has 
a problem at the outermost radii, where the infall 
velocity is much higher than the observed value, 
which is ~ 0.1 km s _1 . One of the main causes for 
the discrepancy is the outer boundary condition 
of the core. In real objects, the infall velocity of 
the outermost region should be zero, which is not 
considered in this asymptotic solution. Hence we 
cannot discuss infall velocities and line profiles in 
our model. Hydrodynamic simulations with re- 
alistic (static) outer boundary conditions show, 
however, that the spatial density distribution and 
temporal variation of the central density of the 
Larson-Penston core are good approximations to 
those in non-magnetized non-rotating cores (Ma- 
sunaga, Miyama, & Inutsuka 1998). 

If a core is supported by rotation, magnetic 
fields, or turbulence, collapse will be delayed, or 
the core will contract slowly. Since molecular 
abundances are determined by a balance between 
chemical and collapse time scales, the abundances 
in such a core are different from those in the 
Larson-Penston solution. In order to investigate 
the effect on molecular abundances, we adopt ana- 
logues of the Larson-Penston core with delayed 
collapse in which the infall velocity is artificially 
decreased by a constant factor /. In comparison 
with observation, we consider the fact that the 
factor / need not be constant but can vary with 
time in reality. In magnetized cores, the collapse 
is slower in earlier, low-density, stages, while in ro- 
tating cores the centrifugal force increases as the 
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collapse proceeds. 

Our chemical reaction network includes gas- 
phase reactions, gas-dust interactions (i.e., adsorp- 
tion and desorption), and diffusive grain-surface 
reactions. These surface reactions occur on grains 
of radius 1.0 x 10~ 5 cm that possess 10 6 bind- 
ing sites. All neutral species are assumed to be 
adsorbed when they hit the grain surface with a 
sticking probability S of unity, unless stated oth- 
erwise. Species on grain surfaces return to the gas 
phase via thermal desorption and impulsive heat- 
ing of the grains by cosmic-rays (Leger, Jura, & 
Omont 1985; Hasegawa & Herbst 1993). For the 
adsorption energies, we adopt those of Model B in 
Paper I; these are basically the values calculated 
by Hasegawa & Herbst (1993) but replaced, if 
available, with experimental data for pure ices. A 
set of grain-surface reactions was taken from Ruf- 
fle & Herbst (2000). To determine reaction rate 
coefficients, those authors generalized the recent 
analysis of Katz et al. (1999), in which a slow 
diffusion rate for H atoms on olivine and amor- 
phous carbon was deduced based on experimental 
measurements of HD formation. Since we prefer 
to consider ice surfaces because these better rep- 
resent the mantles for much of the collapse, we 
calculate reaction rate coefficients based on the 
older approach of Hasegawa & Herbst (1993), in 
which tunneling is included and activation ener- 
gies against diffusion are obtained as a constant 
fraction of adsorption energies. Each surface rate 
coefficient is proportional to the sum of the rates 
for the two reactants to diffuse over the entire 
grain (Hasegawa & Herbst 1993). These rates are 
then altered if necessary according to the revised 
"modified" rate approach of Stantcheva, Caselli, 
& Herbst (2001), which can be described as fol- 
lows. Accretion, evaporation, and diffusion rates 
are calculated for each reactant. If the diffusion 
rate is the smallest among the three values, the 
reaction rate coefficient is unaltered; otherwise the 
diffusion rate in the rate coefficient is replaced by 
the larger of the other two values. The dominant 
alteration is to slow down the reaction rates of 
atomic hydrogen and deuterium. The modifica- 
tion recently proposed for reactions with activa- 
tion energy (Caselli et al. 2002a) was not incor- 
porated. The modified rate approach is a simple 
manner of accounting for the discrete aspects of 
grain chemistry, which are accounted for in greater 



detail by stochastic methods (Stantcheva, Shema- 
tovich, & Herbst 2002). These latter methods are 
not yet available for large reaction networks. 

The gas-grain reaction network has been ex- 
tended to include singly deuterated species. As 
discussed in more detail in §3.7, deuterium frac- 
tionation occurs in both the gas and surface 
phases. Deuterium exchange reactions in the gas 
phase with activation energies are adopted from 
Millar, Bennett, & Herbst (1989). We also tried 
the newly measured and inferred rate coefficients 
for selected fractionation reactions of Gerlich, 
Herbst, & Roueff (2002). The adsorption energy 
and activation barrier against diffusion for surface 
D atoms are assumed to be slightly higher than 
the values for H atoms (Caselli et al. 2002a). In 
total our chemical network consists of 878 species 
and 11774 reactions. In order to check the effect of 
uncertainties in the network, we have also studied 
the core chemistry with the UMIST gas-phase net- 
work (Lc Tcuff, Millar, & Markwick 2000) (§3.6), 
which we have extended to include grain-surface 
reactions and deuterium species in a similar way 
as described above. In this case, 733 species and 
11231 reactions have been considered. 

As in Paper I, we have utilized the so-called 
"low-metal" values for initial gas-phase elemen- 
tal abundances (see Table 1 of Paper I) and have 
adopted the "standard" value of ( = 1.3 x 10~ 17 
s _1 for the ionization rate by cosmic rays. Con- 
sidering that the core is embedded in a molecular 
cloud, A v = 3 mag is assumed at the outer bound- 
ary of the core, so that photodissociation does not 
much affect our results. All heavy elements are as- 
sumed to be initially in atomic form, with carbon 
ionized and oxygen neutral. The initial form of hy- 
drogen is molecular, and the deuterium is assumed 
to be in the form of HD, at a ratio of 3.0 x 10~ 5 
with respect to H2. 

3. Results 

3.1. Molecular Evolution in an Infalling 
Fluid Element 

Figure 1 (a) shows the molecular evolution in 
terms of fractional abundances (n(i)/nn) that oc- 
curs in a fluid element that migrates from a radius 
of 8.2 x 10 3 to one of 1.0 x 10 3 AU in 2.00 x 10 5 yr, 
while the gas density «h increases from 1.7 x 10 4 
to 5.8 x 10 6 cm -3 . Species on the surfaces of dust 
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grains are designated as ices. Only the later stages 
of collapse are depicted because both the density 
variation and molecular evolution are rapid after 
10 5 yr. At this latest stage, the abundances of 
CCS and CO decrease with time while that of 
N 2 H + increases. The reasons for these trends are 
partially the same as described in previous papers 
and partially new; CCS, a so-called "early-time" 
species is depleted via gas-phase reactions and, 
following adsorption, reactions on grain surfaces, 
while CO is adsorbed onto grains, where it reacts 
to produce more complex organic species. The ion 
N 2 H + increases in abundance for at least two rea- 
sons: (i) the binding energy of its parent molecule 
N 2 on grain surfaces is low, which allows its re- 
plenishment into the gas phase, and (ii) its rate 
of depletion is slowed when reactants such as CO 
stick to the grains (Bergin & Langer 1997; Paper 
I). Meanwhile, the abundances of the major ices 
- water, ammonia, and CO - increase with time. 
Methane (CH 4 ) and CO2 ices are formed as well, 
and reach abundances as high as n(i)/nn ~ 10 -5 , 
although they are not shown in the figure. 

Figure 1 (b) shows the molecular evolution oc- 
curring in the same fluid element as in (a), but 
with grain-surface reactions inactivated except for 
H2 formation and ion-electron recombination. The 
results presented are essentially those of Paper I. 
Several major differences in molecular evolution 
are immediately apparent between the result pre- 
sented in Figure 1 (b) and those of our standard 
model (a). First, the ices that are formed mainly 
via surface reactions have greatly reduced abun- 
dances; among these are water, ammonia, carbon 
dioxide and methane. Secondly, carbon monox- 
ide ice, formed mainly via gas-phase reactions fol- 
lowed by accretion, is more abundant in Figure 1 
(b), because the transformation to other species 
via surface reactions is turned off. Thirdly, gas- 
phase species can also have altered abundances. 
The nitrogen molecule has a lowered abundance 
in Figure 1 (b) because its formation on grain sur- 
faces followed by desorption into the gas is signif- 
icant in the normal model (a). The lowered abun- 
dance of molecular nitrogen leads to the low abun- 
dance of N 2 H+. The depletion of CCS is slightly 
accelerated because more oxygen is in atomic form 
rather than in the form of water ice; atomic oxygen 
aids in the transformation of "early-time" species 
such as CCS to CO. 



3.2. Distribution of Molecular Abundances 

The radial distributions of molecular fractional 
abundances within the collapsing core at assorted 
times are obtained by solving for the molecular 
evolution in a similar manner to that described in 
the previous section, but for many fluid elements 
in the flow. The left panels in Figure 2 show dis- 
tributions of gaseous species at t = 1.52 x 10 5 yr, 
1.89 x 10 5 yr, and 2.00 x 10 5 yr (from top to bot- 
tom). At these times, the central density of the 
core is 3 x 10 5 cm -3 , 3 x 10 6 cm~ 3 , and 3 x 10 7 
cm~ 3 , respectively. The distributions of gaseous 
species are qualitatively similar to those of Paper 
I: they are relatively flat at the early stage, labeled 
(a), but the radial dependence becomes more ap- 
parent at later stages, labeled (b) and (c). When 
the central density is 3 x 10 6 cm -3 , NH 3 and N 2 H + 
become more abundant in the inner regions, while 
the other gaseous species are depleted by orders of 
magnitude at the center. When the central density 
reaches 3 x 10 7 cm -3 , even ammonia and N 2 H + 
show a central depletion. 

The right panels in Figure 2 show distribu- 
tions of adsorbed species on the grains. Variations 
among the three evolutionary stages are less pro- 
nounced compared with the gaseous species. Wa- 
ter ice is always the dominant component in the 
grain mantle, while the NH 3 and CH4 ice abun- 
dances are 10-30% relative to water. These three 
species are formed predominantly by surface reac- 
tions in which atomic hydrogen converts adsorbed 
oxygen, nitrogen, and carbon atoms into the sat- 
urated species. Adsorbed CO, formed mainly in 
the gas phase, is transformed to C0 2 , H 2 CO, and 
CH3OH. 

Assuming that the core is spherical, we in- 
tegrate the molecular abundances along the line 
of sight to obtain molecular column densities as 
a function of impact parameter (radial distance) 
from the core center at the three representative 
evolutionary stages. The results are shown in Fig- 
ure 3. While the central density increases from 
3 x 10 5 cm~ 3 to 3 x 10 7 cm~ 3 , the total (hydro- 
gen) column density is enhanced by about an or- 
der of magnitude. The central molecular column 
densities shown in Figure 3, however, show only 
small changes, except for N 2 H + and NH3. The 
column densities of most neutral species, such as 
CO, CS, SO, CCS, and C 3 H 2 , are initially cen- 
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trally peaked, but develop a central hole at later 
stages, while N2H + is centrally peaked even in the 
latest stage, and ammonia shows a very tiny cen- 
tral hole at the latest stage. These qualitative 
features are the same as those reported in Paper 
I. But the column densities of N 2 H+, NH 3 , and 
C3H2 show significant enhancements from those 
in Paper I, which were computed without grain- 
surface reactions. If we consider the time when 
the central density of the core is 3 x 10 6 cm -3 , the 
peak column densities of N 2 H + , NH 3 , and C3H2 
in Paper I are 2.5 x 10 11 cm~ 2 , 3 x 10 13 cm~ 2 , and 
2 x 10 13 cm -2 , respectively, about an order of mag- 
nitude or more lower than reported here. The in- 
crease in the ammonia column density stems from 
the grain-surface formation of N 2 which, unlike 
ammonia ice, is followed by desorption and a se- 
quence of gas-phase reactions leading to gaseous 
NH3. The analogous increase in N 2 H + also fol- 
lows from the grain-surface formation and subse- 
quent desorption of N 2 . Unlike these two cases, 
the enhancement in the column density of C3H 2 , 
an "early-time" species, is the result of a depletion 
of O atoms (§3.1). 

Now that we have summarized the results of our 
standard model, it is useful to see how these results 
depend on some important parameters, such as 
the time scale for collapse, the sticking probability 
on grains, the initial conditions, and the reaction 
network. 

3.3. Collapse Time Scale 

The thick solid lines in Figure 4 show molecular 
column densities as a function of impact parame- 
ter for the case of the Larson-Penston core (/ = 1), 
while cores which collapse more slowly than the 
Larson-Penston model by a factor / of 3 and 10 
are depicted by thick dashed and dotted lines, re- 
spectively. The central densities of the cores are all 
at a density of 3 x 10 6 cm" 3 , so they would appear 
the same if observed via the dust continuum. As in 
Paper I, the molecular column densities are gener- 
ally smaller in the more slowly collapsing cores be- 
cause of heavier depletion (e.g., CCS, CO, SO, and 
CS). The neutral molecule NH 3 and the molecular 
ion N 2 H + show a more complex pattern: here the 
more slowly collapsing core leads to larger column 
densities at larger radial distances but smaller col- 
umn densities near the center of the core. At outer 
radii, the parent molecule N 2 is more abundant in 



the slow collapse models because of its slow forma- 
tion rate, while at inner radii, on the other hand, 
more nitrogen is depleted onto grains as NH 3 ice 
in the slow collapse models. 

In general, the sensitivity of the column densi- 
ties to the collapse time scale is smaller than that 
obtained in Paper I. For example, the central col- 
umn densities of CCS and C3H 2 in the / = 3 
model are lower than those in the / = 1 model 
by more than an order of magnitude in Paper I, 
whereas here they are lower by factors of 6 and 2, 
respectively. What causes this change? The abun- 
dances of the early-time species, which are typi- 
cally unsaturated (H-poor) organic molecules such 
as CCS and C3H2, are sensitively dependent on 
the C/O elemental ratio in the gas phase (Pratap 
et al. 1997; Terzieva & Herbst 1998) according 
to gas-phase models of cold clouds because this 
elemental ratio determines how much atomic oxy- 
gen is available in the gas. In the current models, 
early-time species are relatively abundant even in 
the slowly collapsing cores because of the depletion 
of O atoms via water formation on grain surfaces. 
The dependence of the NH3 and N 2 H+ column 
densities on the collapse time scale is also smaller 
than in Paper I, where the central densities of NH 3 
and N 2 H+ in the / = 10 model are larger than 
those in the / = 1 model by factors of 5 and 10, 
respectively, because of their slow formation in the 
gas phase. Here, the major synthetic pathways for 
both species lie through N 2 that is desorbed from 
grain surfaces. The increase in time available due 
to slow collapse has less of an effect. 

3.4. Sticking Probability 

The thin lines in Figure 4 show molecular col- 
umn densities vs. radial distance for slowly col- 
lapsing cores with low sticking probabilities S. 
Specifically, for an / of 3 (dashed lines), S is 0.33, 
and for an / of 10 (dotted lines), S is 0.1. The 
lowered sticking probabilities are chosen to com- 
pensate for the increased time spent in collapse. 
The compensation is of course only a partial one, 
because rate equations are non-linear, including 
many reactions other than adsorption. 

The dependence on the sticking probability 
varies among species. The column density of CO 
is higher in models with a lower sticking proba- 
bility. Grain-surface reactions transform CO to 
C0 2 , H 2 CO, and CH 3 OH, which mostly remain 
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on the grain because of higher adsorption ener- 
gies. In the models with lower sticking probabil- 
ity, such processes are less efficient. The column 
density of SO is higher in models with a lower 
sticking probability as well, because depiction of 
SO itself is slow enough to allow the late peak to 
grow after the depiction of carbon atoms, which 
are the dominant reactant during the early-time 
chemistry. Although the gas-phase abundances of 
CCS, CS, and C3H2 are highest in the inner re- 
gions for the standard Larson-Penston case, for 
models with / > 1 the abundances are higher in 
models with a lower sticking probability and thus 
with a smaller rate of depletion. At larger radii, 
on the other hand, C3H2 is less abundant in the 
models with a lower sticking probability, because 
of the higher abundance of destructive O atoms 
in the gas phase. Gas-phase N 2 and NH 3 tend 
to increase in abundance at small radii with lower 
sticking rates, because the lessened rate of deple- 
tion via grain adsorption is more important than 
the lessened rate of formation, while N 2 H + is less 
abundant because of the higher abundance of CO, 
reactions with which deplete the ion. 

3.5. Initial Conditions 

So far we have assumed that the initial central 
density of the core, nn(t — 0, R = 0), is 2 x 10 4 
cm" 3 , which is often smaller than observed in qui- 
escent dense cores but is often found in molecular 
clouds (e.g. Onishi et al. 1996). Since the initial 
condition of the core is uncertain, it is interest- 
ing to determine the dependence of the molecular 
column densities on the initial gas density. Fig- 
ure 5 shows molecular column densities in cores 
with initial central densities of 2 x 10 4 (solid lines), 
2 x 10 3 (dashed lines), and 2 x 10 2 cm" 3 (dotted 
lines). In all three models, the initial gas consists 
of our standard conditions: atoms and ions except 
for hydrogen, which is in molecular form. The cen- 
tral density in all models is 3 x 10 6 cm" 3 , which 
is reached after 1.89 x 10 5 yr, 6.33 x 10 5 yr, and 
2.04 x 10 6 yr, respectively. It should be noted that 
the latter two models are in a low density state 
for a significant fraction of the time; their central 
densities reach 2 x 10 4 cm" 3 at t = 4.44 x 10 5 
yr and 1.85 x 10 6 yr, after which the evolution of 
the density distribution is the same as that of the 
first model with the time t of calculation shifted 
accordingly. 



It can be seen in Figure 5 that the dependence 
on the initial conditions is a factor of ~ 2 for most 
of the species, which suggests that the early evo- 
lutionary stages with low densities do not much 
affect the molecular abundances in the dense core. 
Exceptions are N 2 H + and NH 3 , which are less 
abundant in the models with a low initial den- 
sity. In these models, the higher abundance of H 
atoms enhances the surface formation of NH3 at 
the expense of N 2 , which is the major precursor 
of gaseous ammonia as well as N 2 H + , since sur- 
face ammonia does not desorb efficiently in our 
model. In reality, the dependence of the N 2 H + 
and NH3 abundances on the initial density could 
be smaller. In clouds with low density (~ 10 2 
cm" 3 ) the temperature would be higher than as- 
sumed here (10 K), and hence efficient desorption 
of H atoms would suppress the surface formation 

of NH3. 

3.6. Uncertainties in Reaction Networks 

Molecular evolution proceeds via a chemical re- 
action network that consists of a large number of 
reactions both in the gas phase and on the grain 
surfaces. The gas-phase chemistry has been stud- 
ied by more investigators over a longer period of 
time, and is generally better understood than the 
grain-surface chemistry. Gas-phase models have 
succeeded in reproducing or predicting the abun- 
dances of most (80 %) of the observed species 
in quiescent sources, in understanding four-order- 
of-magnitude deuterium enrichments, and in ex- 
plaining the existence of peculiar species such as 
carbon-chain molecules (Terzieva & Herbst 1998; 
Millar, Bennett, & Herbst 1989; Roberts & Mil- 
lar 2000). Yet, in spite of these successes and 
significant efforts in refining the networks, un- 
certainties in the chemistry still remain. These 
uncertainties arise from a variety of causes. First, 
most reactions in the models have not been stud- 
ied in the laboratory and their rates are esti- 
mated based on analogous better-studied sys- 
tems. These estimates are on somewhat firmer 
ground for ion-molecule reactions than they are 
for neutral- neutral reactions, since it is always 
difficult to determine whether the latter possess 
significant activation energies. The extent of the 
role of neutral-neutral reactions remains a signif- 
icant problem (Terzieva & Herbst 1998). Sec- 
ondly, most applications of the chemical models 
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are at low temperatures (e.g., 10 K) whereas most 
laboratory experiments are at higher tempera- 
tures. The temperature dependence can often be 
obtained from theoretical considerations or com- 
parison with other systems measured over a large 
temperature range, but not always with certainty. 
A simple example concerns the rate coefficients 
for ion-polar neutral reactions, for which some au- 
thors include temperature- independent rate coeffi- 
cients while others, extrapolating from theory and 
a few measured results, include an inverse tem- 
perature dependence. The strong but uncertain 
temperature dependence of radiative association 
reactions also leads to significant uncertainties. 
Thirdly, even subtle effects such as the role of fine 
structure can make a conversion from a laboratory 
rate to an interstellar rate rather uncertain. 

Currently there are two partially-independent 
extensive gas-phase chemical reaction networks: 
the UMIST data base (Le Teuff, Millar, & Mark- 
wick 2000) and the New Standard Model (NSM) 
(Terzieva & Herbst 1998). Since the authors of 
these networks have different policies and stan- 
dards in compiling the data, the differences in 
the two models can give us a sense of how uncer- 
tainties in chemistry affect the results. But this 
is only a sense, because the absolute uncertain- 
ties in reaction rate coefficients probably exceed 
in many instances the differences between the two 
networks. Markwick (2002) has used a Monte 
Carlo procedure to show how the estimated uncer- 
tainties in the latest rates adopted by the UMIST 
compilers lead to uncertainties in the predicted 
abundances, an effect which becomes significant 
for large molecules. 

We have adopted the NSM for the gas-phase 
reactions in our calculation. Here we replace it 
with the UMIST data base, and check how sig- 
nificantly our results are affected by uncertain 
rate coefficients. The dash-dotted lines in Figure 
5 show the molecular column densities obtained 
with the UMIST data base, which we extended 
to include the same grain-surface reactions and 
gas-grain interactions as in our NSM-based model. 
The physical model is the Larson-Penston flow 
with an initial central density of 2 x 10 4 cm~ 3 . 
Comparison with results obtained with the NSM 
(solid lines in Figure 5) suggests that the uncer- 
tainties in the gas-phase network cause errors of 
a factor of < 2 in the column densities of most 



molecular species. Exceptions are sulphur-bearing 
molecules; the CCS and CS column densities are 
higher with the UMIST data base than with the 
NSM by factors of 4 and 30, respectively. These 
differences arc not surprising because it has been 
known for some time that the differences among 
S-bearing species exist even in the purely gas- 
phase pseudo-time dependent (i.e. constant den- 
sity) models. In an effort to locate the specific 
source of the discrepancies, we checked the main 
formation and destruction paths for each molecu- 
lar species and found several important reactions 
which have different rate coefficients in the two 
networks. 

One of the critical reactions is O + HCS + , since 
this reaction depletes protonated CS and inter- 
feres with the production of CS via dissociative 
recombination with electrons. The rate coefficient 
for this reaction is 5.0 x 10~ 12 cm 3 s _1 in the 
UMIST data base and 5.0 x lO" 10 cm 3 s" 1 in the 
NSM. It is of interest to illuminate how this dif- 
ference came about. Since the UMIST website 
(http://www.rate99.co.uk/) shows how the rate 
coefficients for individual reactions are obtained, 
one can see that the rate coefficient for O + HCS + 
stems from an early paper by Millar & Herbst 
(1990) in which it is actually written that "both 
rate coefficients and products are uncertain." The 
NSM value is presumably based on other stud- 
ies involving atomic oxygen, which tend to show 
rate coefficients up to one order-of-magnitude be- 
low the Langevin rate. Another critical reaction is 
that between C and SO. The UMIST value stems 
from a simple theoretical model used in the shock 
study of Lccn & Graff (1988) while the origin 
of the much lower NSM value is uncertain. From 
these examples, one can see that different choices 
can be made when modelers are faced with less 
than complete information. Regarding the very 
uncertain rates of neutral-neutral reactions, one of 
us (E.H.) has recently started to collaborate with 
a European network of chemical kineticists to ob- 
tain expert opinion on unstudied systems. The 
resulting rate coefficients for neutral-neutral reac- 
tions will no longer contain what once was the 
standard temperature dependence of (T/300) 0,5 
derived from hard sphere arguments. 
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3.7. D/H ratio 

The elemental abundance of deuterium is ~ 
10" 5 relative to hydrogen in the local interstellar 
medium. In molecular clouds, deuterium is pre- 
dominantly in the form of HD, but is transferred 
from this reservoir to other molecular species by 
proton/deuteron exchange reactions, the most im- 
portant of which is H+ + HD -> H 2 D+ + H 2 
(Millar, Bennett, & Herbst 1989; Roberts & Mil- 
lar 2000). The exchange reactions are exother- 
mic by amounts of energy large enough that the 
backwards endothermic reactions are inefficient in 
clouds with temperatures significantly lower than 
the reaction energy expressed as a temperature. 
As a result, D/H ratios such as H 2 D + /H^ are 
much larger than the elemental D/H ratio, an 
effect known as fractionation. The deuterated 
isotopomers formed in exchange reactions react 
with other species to transfer the deuterium fur- 
ther, producing a significant deuterium fractiona- 
tion for many minor species in the model. The ion 
H 2 D + is depleted at low temperatures principally 
by reactions with electrons and neutral species 
with high proton affinity such as CO and H 2 0. 
If these heavy species are depleted from the gas 
onto the dust grains, then the D/H ratios in H^ 
and in other molecular species are enhanced if, as 
is the case here, the ionization degree is low (i.e. 
n(e)/riH < 3 x 10" 7 ). These arguments suggest 
that D /H ratios of species other than HD can be 
used as probes of cloud conditions such as temper- 
ature, degree of ionization, and degree of molec- 
ular depiction. It is therefore useful to present 
molecular D/H ratios for our models. 

Fractionation can also occur via grain-surface 
chemistry. In particular, D atoms are formed in 
the gas phase via the dissociative recombination 
of molecular ions such as H 2 D + and DCO + , and 
then are adsorbed onto grain surfaces, where they 
can react in much the same manner as atomic hy- 
drogen. For example, if a surface CO molecule is 
subjected to reactions with atomic H and D, both 
the normal species H 2 CO and CH 3 OH can be pro- 
duced as well as their singly and multiply deuter- 
ated isotopomers (Charnley et al. 1997; Caselli et 
al. 2002a). 

Figure 6 shows radial distributions (n(i)/nn) 
of deuterated and normal species in the Larson- 
Penston core when the central density of the core 



is 3 x 10 5 cm" 3 (top), 3 x 10 6 cm" 3 (middle), 
and 3 x 10 7 cm" 3 (bottom). The left panels 
show gas-phase species, and the right panels show 
abundances in the grain mantles. As expected, 
the D/H ratios in the gas-phase molecules are 
higher at smaller radii, and at later evolution- 
ary stages because of the stronger depletion of 
heavy molecules and the lowered ionization de- 
gree. The ionization degree obtained in the cal- 
culation can be approximately fitted by the sim- 
ple function 6 x 10 _9 (10 5 cm _3 /riH) ' 5 in regions 
with relatively high density (rin > 10 5 cm" 3 ) for 
the Larson-Penston model, and varies by a fac- 
tor of < 2 depending on the parameter / and S. 
This degree of ionization is close to what one esti- 
mates if all ionization stems from the interaction 
of cosmic rays and molecular hydrogen, and is also 
consistent with the value (~2x 10" 9 ) estimated 
by Caselli et al. (2002b) at the center of L1544 
where nn ~ 10 6 cm" 3 . The resulting deuterium 
fractionation can be so high that the deuterated 
isotopomers can become more abundant than the 
normal isotopomers; H 2 D + is more abundant than 
H;j~ , and is the most abundant ion in the innermost 
regions when the density is high. Given the high 
degree of deuterium fractionation, further studies 
with multiply deuterated isotopomers would be in- 
teresting in relation to the recent observation of 
these species in protostars. 

Among the surface species, the radial distribu- 
tions show somewhat less of a dependence on den- 
sity and radius. Still, there are exceptions, par- 
ticularly singly deuterated formaldehyde - HDCO 
- which, for the higher density inner cores, ac- 
tually becomes more abundant than the normal 
isotopomer. On the other hand, the more abun- 
dant of the two singly deuterated isotopomers of 
methanol - CH 2 DOH - only achieves an abundance 
approximately 0.01 of normal methanol. 

So far, we have discussed the radial distribu- 
tion in terms of fractional abundances and not 
the column densities of deuterated and normal iso- 
topomers as functions of impact parameter. Fig- 
ure 7 shows column density ratios for some deuter- 
ated species to normal species as a function of im- 
pact parameter from the core center when the cen- 
tral density is 3 x 10 6 cm" 3 . The Larson-Penston 
core is represented by solid thick lines, while col- 
lapses slowed down by factors / of 3 and 10 are 
depicted by thick dashed and dotted lines, respec- 
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tively. For these latter two cases, models were also 
run with sticking probabilities of 1//; the results 
are, as usual, shown as thin dashed and dotted 
lines. As can be seen, the D/H ratios are higher in 
models with heavier depiction; i.e. those with slow 
collapse and a high sticking probability. The ra- 
tio, which is an average along the line of sight, also 
varies with species; it is higher if the radial distri- 
bution of the species is more centrally peaked be- 
cause deuterium fractionation is greater in the in- 
ner regions. For example, the molecular ion N 2 H + 
and the "late-time" species NH 3 have higher ratios 
than neutral "early-time" species such as HNC. 

In the models presented above, the rate coef- 
ficients for the deuterium fractionation exchange 
reactions are adopted from Millar, Bennett, & 
Herbst (1989). Recently, Gerlich, Herbst, & 
Roueff (2002) and Gerlich & Schlemmer (2002) 
measured these rate coefficients in an ion trap at 
10 K, and found that the rate coefficients of the 
forward reactions (e.g. fht + HD -> H 2 D+ + H 2 ) 
are smaller and those of backward reactions prob- 
ably larger than those utilized by Millar, Bennett, 
& Herbst (1989). Roberts et al. (2002) showed 
that these revisions of the rate coefficients signifi- 
cantly lower the D/H ratios in quiescent cores for 
both gas-phase and accretion-containing models. 
To determine the effect of the revised rates for 
prestellar cores, we used theses rates (summarized 
in Tabic 1 of Roberts et al. (2002)) with a core 
undergoing Larson-Penston collapse. Some results 
are shown with dashed-dotted lines in Figure 7. 
The D/H ratios are typically lowered by an order 
of magnitude. 

4. Discussion 

4.1. Comparison with L1544 and Previous 
Models 

The prestellar core LI 544 is the most intensely 
studied object of this class. Its radius is about 
15000 AU, and the central density n(H 2 ) is esti- 
mated to be 1.5 x 10 6 cm~ 3 from dust continuum 
measurements. Various molecular lines have been 
observed and gas-phase molecular column densi- 
ties estimated. Detailed mapping and interferom- 
eter observations have revealed that CO and CCS 
are depleted at the core center, while the N 2 H + 
emission peaks at the center (Caselli et al. 1999; 
Ohashi et al. 1999; Tafalla et al. 1998; Williams 



et al. 1999). 

Firstly, we compare molecular column densi- 
ties obtained in our models with those observed at 
LI 544. In Table 1, we list the estimated molecu- 
lar column densities from observational data which 
are obtained at the core center. Table 1 also lists 
molecular column densities calculated in our mod- 
els with a variety of / and S values when the cen- 
tral density is hh = 3 x 10 6 cm~ 3 . The tabu- 
lated results are surface-averaged (i.e. weighted 
by 2irRdR) within a radius of 1000-2000 AU de- 
pending on the actual beam size of each observa- 
tion. The bold numbers indicate agreement with 
the observation to within a factor 3, while the 
numbers in italics indicate disagreement by a fac- 
tor of > 5. The model with / = 1 shows the 
best agreement, while the very slow collapse model 
(/ = 10, S = 0.1) shows the worst agreement. 

The factor /, the ratio of the actual collapse 
time scale to the Larson-Penston collapse time 
scale, should be considered with caution. Al- 
though we set / to be constant for simplicity, in re- 
ality it varies with time. If the core is supported by 
magnetic fields, / is large in the early low-density 
stages, but decreases to unity when the core be- 
comes too dense to couple with magnetic fields. 
On the other hand, / increases in the later stages 
of collapse via rotational support, if the core has 
significant angular momentum. Since chemical re- 
actions generally proceed faster at higher densi- 
ties, the value of the factor / at higher densities 
is more important in determining the molecular 
abundances in cores. Therefore the fact that the 
worst agreement is obtained with the highest value 
of / indicates that the collapse time scale at a cen- 
tral density of ~ 10 6 cm~ 3 should not be larger 
than the free-fall time scale by a factor of 10. 

Secondly, we compare the radial distribution 
of fractional molecular abundances in our model 
cores with those in L1544. Tafalla et al. (2002) 
utilized radiative transfer calculations to estimate 
the distribution of fractional molecular abun- 
dances in prestellar cores such as L1544 from 
observational data. These distributions, which 
represent a new development, arc not to be con- 
fused with column density measurements. Figure 
8 (a), taken from formulae of Tafalla et al. (2002), 
shows their radial distribution of molecular abun- 
dances in L1544. At the core center, CO and 
CS are heavily depleted and NH 3 is slightly cn- 
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hanced, while the abundance of N 2 FI + is constant. 
Figures 8 (b) and (c) show our results for the L-P 
core (/ = 1,5 = 1) and the model with / = 3 
and 5=1 when the central density of the core 
is 3 x 10 6 cm~ 3 . The calculated distributions are 
in qualitative agreement with those determined 
in L1544, but not in quantitative agreement, ex- 
cept for NH 3 . In our models, the N 2 H+ fractional 
abundance is significantly higher towards smaller 
radii, while it is inferred to be constant in L1544. 
Also, the gradient of the CO abundance in the cen- 
tral region of LI 544 is steeper than in our models. 
Although the steep CO gradient in L1544 is partly 
caused by the exponential function assumed in the 
analysis, the CO intensity profile obtained from 
our L-P core is still more centrally peaked than in 
L1544 (Tafalla 2003, private communication). In 
the slow collapse model (c), CO is depleted over a 
larger area, so that the gradient of the CO abun- 
dance is even shallower than in the L-P model. 

The centrally peaked N 2 H+ abundance distri- 
bution in our models is caused by the relatively 
slow formation of N 2 in the outer regions, which 
are at low density. The degree of CO depletion, 
on the other hand, is dependent on its adsorption 
energy (Paper I). To determine if the agreement 
with the results of Tafalla et al. (2002) can be im- 
proved, we performed new calculations in which 
the nitrogen is all in the form of N 2 initially and 
the adsorption energy of CO on grains is varied. 
In our models discussed up to now, the adsorp- 
tion energy for CO is that of pure CO ice, which 
is 960 K. The CO adsorption energy is higher on 
other surfaces; e.g., it is 1210 K if CO is adsorbed 
on silicate material and 1780 K on water (polar) 
ice (Paper I and references therein). We have run 
models with both 1210 K and 1780 K for the ad- 
sorption energy of CO. The resulting distributions 
are shown in Figure 8 (d-f). As can be seen, the 
steep drop in CO abundance is well reproduced 
in the model with / = 3 and an adsorption en- 
ergy of 1780 K (f). Concerning the distribution of 
N 2 n + , the constant abundance observed in LI 544 
is better reproduced with the N 2 -rich initial abun- 
dance, although the initial fraction of N 2 should 
be smaller than unity in order to reproduce the ab- 
solute abundance. The improvement in the calcu- 
lated N 2 H + distribution suggests either that early 
formation of N 2 is more efficient than assumed in 
our standard model, or that some molecular devel- 



opment occurs before the collapse stage of prestcl- 
lar cores. 

How does our current agreement with observa- 
tion in L1544 compare with those based on previ- 
ous models ? In Paper I we preferred the Larson- 
Penston model over models with slower collapse. 
With this model, we were able to best reproduce 
the radial distributions of the column densities of 
CCS, CO, and N 2 H+ as well as the peak col- 
umn densities of CCS and CO. Our calculated 
column density for N 2 H+ was low by a factor of 
25, however. The current results are a consider- 
able improvement since they show agreement over 
a wider group of species, which includes N 2 H + . 
This agreement is obtained without slowing down 
the core collapse. Moreover, we see no improve- 
ment when the collapse is retarded except for fit- 
ting the results of Tafalla et al. (2002), which also 
requires some fraction of nitrogen to be initially in 
its molecular form and CO to have a much higher 
adsorption energy. The improvement from Paper 
I is thus due to the inclusion of grain chemistry. 

Li et al. (2002) include surface chemistry in 
their model. Their chemistry is not identical with 
ours because they use a small subset of the UMIST 
gas-phase network and our grain surface reactions. 
With their networks and the absence of a mag- 
netic field, they obtain reasonable agreement with 
the data except that CS and CCS are found to 
be centrally peaked, in disagreement with obser- 
vation and also in disagreement with our results. 
Considering the results of §3.6, the difference be- 
tween our models could be caused by the use of 
different chemical networks. When a magnetically 
retarded collapse is used, they find a spatial hole at 
the core. Their magnetic model is best compared 
with our model in which / — 3 and 5 = 0.33 since 
in the work of Li et al. (2002) the sticking proba- 
bility is assumed to be 0.3, and the collapse time 
scale of the magnetized cloud is larger than that 
of non-magnetized cloud by a factor of ~ 2 when 
the central density is 10 5 — 10 6 cm~ 3 . Agreement 
between our model with 5 = 3 and / = 0.33 and 
observation is improved from Paper I by the inclu- 
sion of the grain-surface reactions, although it is 
certainly not our best model. Although some dif- 
ferences still remain between this current model 
and Li et al. (2002), we can conclude that the 
main cause of the difference between our Paper I 
and Li et al. (2002) is grain-surface chemistry. 
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Table 1 

Molecular Column Densities (cm -2 ) in L1544 and Theoretical Models 



Species 


LI 544 


f=l, S=1.0 


f=3, S=1.0 


f=3, S=0.33 


f=10, S=1.0 


f=10, S=0.1 


CO 


1.2(18) a 


1.1(18) 


6.1(17) 


2.0(18) 


3.3(17) 


3.0(18) 


CCS 


2.6(13) b 


1.9(13) 


3.0(12) 


9.6(12) 


2.0(11) 


5.7(12) 


N 2 H+ 


2.0(13) a 


1.3(13) 


7.7(12) 


2.7(12) 


7.5(12) 


3.3(12) 


N 2 D+ 


4.3(12) a 


5.1(12) 


2.7(12) 


6.6(11) 


2.0(12) 


4.8(11) 


NH 3 


4.0(14)° 


3.6(14) 


2.3(14) 


2.2(14) 


2.0(14) 


3.6(14) 


CS 


4.6(13) c 


4.7(13) 


1.5(13) 


4.1(13) 


2.3(12) 


2.8(13) 


HNC 


1.5(14) d 


3.7(14) 


1.6(14) 


1.7(14) 


8.6(13) 


2.8(14) 


HCO+ 


1.1(14)" 


1.3(14) 


1.3(14) 


1.8(14) 


1.0(14) 


1.7(14) 


DCO+ 


4.0(12) a 


1.8(13) 


2.6(13) 


2.8(13) 


1.7(13) 


3.2(13) 



a Caselli et al. (2002b) 
b Ohashi et al. (1999) 

c Calculated from best fit model of Tafalla et al. (2002), assuming a spherical core with radius of 
15000 AU. 

d Hirota et al. (2003) 



4.2. Variation among Cores 

Several groups have performed surveys of dense 
cores through various molecular lines (Suzuki et 
al. 1992; Hirota et al. 1998; Tafalla et al. 2002; 
Caselli et al. 2002b). Comparison of the results 
of our model with these surveys should be use- 
ful in investigating how molecular column densi- 
ties vary during the evolution of cores. Table 2 
lists observed column densities of assorted species 
in 10 prestellar cores. We chose compact (~ 10 4 
AU) prestellar cores, which, we believe, have just 
begun collapse or are on the verge of star forma- 
tion. We carefully selected the data obtained to- 
wards the peak position of N 2 H + emission, in or- 
der to avoid contamination with the spatial vari- 
ation within the cores. The column densities in 
Table 2 with parentheses arc observed somewhat 
(<~ 30") offset from the N 2 H + peak. Since different 
beam sizes can also cause variations in estimated 
column densities, we chose observations with sim- 
ilar beam sizes for each molecular species. If the 
core has been observed and analyzed by Tafalla et 
al. (2002), we calculate from their best-fit model 
the surface-averaged column density within the 
beam size of Caselli et al. (2002c) (54") for N 2 H+, 
Suzuki et al. (1992) (80" and 62") for NH 3 and 
CCS, and Hirota et al. (1998) (34") for CS. The 
column density of CO is averaged over the central 
radius of 2000 AU. The molecular column densi- 
ties of L1544 listed in Table 2 are slightly different 



from those in Table 1, because of the larger beam 
sizes adopted in Table 2. Judging from the errors 
explicitly given for N2H + and differences in esti- 
mated column densities by independent observers, 
the listed values may contain uncertainties by a 
factor of <~ 2 — 3. Large variations in the col- 
umn densities of certain species are still apparent 
among cores, especially among the four objects at 
the top and bottom of Table 2. Specifically, in 
L1521E and L1521B, column densities of N 2 H+ 
and NH 3 are smaller and those of CCS and CS 
are larger by about an order of magnitude than in 
L63 and TMC-2A. 

Let us first consider how our results depend 
on the central density (nn) of the core. Figure 
9 shows calculated molecular column densities as 
a function of the central density in our models 
when the central density is nn > 3 x 10 4 cm~ 3 . 
The column densities of CCS, N 2 H+, NH 3 , and CS 
are surface-averaged within the typical beam sizes 
mentioned above at the core center, while those 
of CO and HCO + are surface-averaged within a 
radius of 2000 AU. The column densities of N 2 H+ 
and NH 3 tend to increase with increasing central 
density; i.e., they increase as the core evolves. The 
ion N 2 H + is of special interest because its column 
density increases monotonically in all of our mod- 
els, including the models with an initially high 
abundance of N 2 discussed in §4.1. In addition, 
the N 2 H + column density does not depend signifi- 
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cantly on the total mass (size) of the core, because 
its absolute abundance is higher in the inner re- 
gions. Hence this ion appears to be a useful probe 
of core evolution. On the contrary, CCS and CS 
decrease with time in most of the models, and CO 
and HCO + show little temporal variation except 
for the very early stages, when tih < 10 5 cm~ 3 . A 
variation among core models is apparent for CCS, 
CO, and CS. 

In Table 2 we have arranged the objects in the 
order of increasing N 2 H+ column density, assum- 
ing this order to represent evolutionary stage. We 
see that Table 2 indicates that the column den- 
sity of NH 3 tends to increase with increasing time, 
which is consistent with our model results, as is the 
decrease in CCS and the weaker but similar trend 
for CS (see Figure 9). Another interesting fact is 
that in L1521E and L1521B, which are considered 
to be the youngest ones in Table 2 according to 
the N 2 H + column density, CCS emission peaks at 
the core center (Ohashi 2000; Hirota ct al. 2002). 
Such a CCS distribution is indeed characteristic of 
a very early stage of core evolution in our models. 

Since the variation among core models is rela- 
tively significant for CCS, CO, and CS, a quanti- 
tative comparison with observation might exclude 
some collapse models, assuming that the collapse 
time scale (i.e. the parameter /) is the same in all 
cores. For example, the computed temporal vari- 
ation of CCS is much smaller in the L-P model 
than the observed variation in Table 2. However, 
the morphology of the core should be considered in 
a quantitative comparison of these species, which 
are heavily depleted in the central region. Our 
caution arises because the temporal variation of 
their model column densities is small at least par- 
tially because of the spherical symmetry assumed 
in our calculations. Although their abundances 
decrease significantly in the central regions as the 
core evolves (Figure 2), contributions from the 
outer radii keep the molecular column densities al- 
most constant. In reality, cores are not spherical, 
and depletion in the central region will be more 
apparently reflected in the molecular column den- 
sity. Models of non-spherical cores are desirable 
for further studies. 

LI 544 has significantly higher column densities 
of CCS and CO compared with other objects with 
a similar N 2 H+ column density. There are a few 
possible explanations. Firstly, the difference can 



be caused by core geometry and/or view angles. If 
the cores are not spherical and are observed with a 
nearly face-on angle, we would obtain significantly 
lower column densities of CCS and CO than in 
the case of a spherical core or a non-spherical core 
with edge-on geometry, as discussed above. For 
this explanation to be germane, L1544 would be 
a spherical core or a core with edge-on geometry, 
while the others would be non-spherical cores with 
a more face-on geometry. Interestingly, Ohashi et 
al. (1999) suggested from their CCS observation 
that L1544 is a non-spherical core that is observed 
with a nearly edge-on geometry. Ciolek & Basu 
(2000) expected a similar geometry by comparison 
of the object with their magnetized core model. 
Another possible explanation is that the collapse 
time scale is different among cores. While LI 544 
is best reproduced by the Larson-Penston core in 
our model, other objects may collapse more slowly 
(e.g. by ambipolar diffusion). 

4.3. D/H Ratios In Assorted Cores 

Finally, we compare the molecular D/H ratios 
obtained in our models with observations in a va- 
riety of cores. As mentioned in §3.7, molecular 
D/H ratios are predicted to increase during the 
evolution of the cores. It would be useful to plot 
these ratios vs the central densities, which are a 
direct measure of the evolution of the cores. Since 
we do not know the central density (i.e. evolu- 
tionary stage) of most evolving cores, we can use 
the column density of the ion N 2 H + as an indica- 
tor of evolution, following our previous arguments 
in §4.2. This ion has been studied in approxi- 
mately 60 cores by Caselli et al. (2002c), while 
the column density ratio DCO + /HCO + has been 
measured by Butner et al. (1995) and Williams 
et al. (1998), and the ratio DNC/HNC has been 
measured by Hirota et al. (1998) and Hirota et 
al. (2003), both towards several prestellar cores. 
The ratio of N 2 D+ /N 2 H+ in L1544 has been stud- 
ied by Caselli et al. (2002b). We emphasize that 
the column density of N 2 H+ does not depend sig- 
nificantly on the total mass of the core, because 
of the centrally peaked distribution of its absolute 
abundance. The D/H column density ratios are 
not significantly dependent on the total mass ei- 
ther, to the extent that the dependence is at least 
smaller than that of the molecular column densi- 
ties themselves. Thus, it is appropriate to include 
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Table 2 

Observed Molecular Column Densities in Assorted Prestellar Cores 



Object 


N 2 H+ 
10 12 cm" 2 


NH 3 
10 14 cm- 2 


CCS 
10 12 cm" 2 


CS 
10 14 cm- 2 


CO 
10 17 cm" 2 


L1521E 


< 0.14 a 


0.73 a 


28 a 


3.0 a 




L1521B 


l±0.85 b 


0.6° 


36 c 


1.3 d 




L1517B 


3±0.3? 3.1 f 


7? 2.1 f 


8.6° 


0.13 f 


1.5 f 


L1400K 


4±2? 2.7 f 


1.8? 7.2 f 


3.6 C 


0.25 f 


2.1 f 


L1512 


3±2° 


( 7 )c,g 


(2.9) c 






L1498 


8± 4? 3.0 f 


(4.1)? 2.3 f 


(16.5) c 


(0.40)?0.14 f 


2.4 f 


L1495 


6.0 f 


2.4 f 




0.13 f 


3.6 f 


LI 544 


9±2? 7.3 f 


1.8 f 


20 h 


0.46 f 


4.0? 12 s 


L63 


8±4° 


7.9° 


< 1.7 C 


0.36 d 




TMC-2A 


11±3° 


10.7 C 


3.2 C 







a Hirota ct al. (2002) 
b Takakuwa ct al. in prep 
°Suzuki ct al. (1992) 

d Hirota ct al. (1998), assuming 34 S/S= 4.2% 
"Caselli et al. (2002c) 
f Tafalla ct al. (2002) 

s The column densities with parentheses are observed somewhat (~ 30") offset from 
the N 2 H+ peak. 

h Ohashi ct al. (1999) 

'Casclli et al. (2002b) 
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cores of differing mass on the same evolutionary 
plot. 

Figure 10 contains our plot of D/H column den- 
sity ratios for HNC, HCO + and N2H + versus the 
N 2 H + column density. Again we chose the ob- 
servational data which are measured at or close to 
the emission peak of N 2 H+: L1498, L134A, L1544, 
L1696A, L183(S), L63, and L1155C for the plot 
of DCO+/HCO+, and L1521B, L1082A, L1521E, 
L1512, L1544, L63, and L1155C for DNC/HNC. 
Figure 10 shows that D/H ratios tend to increase 
with N 2 H+ column density, which is in quali- 
tative agreement with the model results. For 
each molecule, the results from three different gas- 
phase chemical networks are depicted in order to 
consider uncertainties arising from uncertainties 
in the chemistry. In addition to the NSM and 
UMIST networks discussed earlier, we utilize the 
NSM network with the new rates for deuterium 
fractionation of Gerlich, Herbst, & Roueff (2002) 
as utilized by Roberts et al. (2002), which we la- 
bel RHM. The theoretical results, plotted as lines, 
encompass our model cores as their central den- 
sity increases from 3 x 10 5 cm~ 3 to 3 x 10 7 cm~ 3 . 
The calculated molecular column densities are av- 
eraged over radii of 1400 AU, 4900 AU, and 3780 
AU from the core centers for HNC (DNC), HCO+ 
(DCO + ), and N 2 H + , respectively, considering the 
beam sizes of the observations. In calculating the 
ratio of N 2 D + /N 2 H + , we averaged the column 
density over a radius of 1400 AU, referring to the 
beam size of Caselli et al. (2002b). For each D/H 
ratio and model network, we calculated results for 
the Larson-Penston case, cases with collapse de- 
layed by factors / of 3 and 10, and cases in which 
a sticking probability of 1/f is assumed. As dis- 
cussed, the theoretical plots tend to show a direct 
relation between the D/H ratios and the N 2 H+ 
column density. The ratio of DNC/HNC is of spe- 
cial interest because its dependence on the collapse 
model is relatively significant. 

If we initially restrict our attention to the NSM 
calculations, we see that the standard Larson- 
Penston solution is the best overall representation 
of all but the largest DNC/HNC value, which oc- 
curs for L63. This agreement, however, is not 
maintained with the RHM network, in which less 
efficient fractionation must be balanced by the 
adoption of very slow collapse models. Even 
so, the RHM model cannot account for the high 



N 2 D+/N 2 H+ ratio observed in L1544. The re- 
sults with the UMIST network also vary somewhat 
from those of the NSM network. In summary, we 
can state that despite all of the uncertainty, the 
use of the new RHM values makes accounting for 
deuterium fractionation in a quantitative manner 
far more difficult, a result previously obtained by 
Roberts et al. (2002) for static cores. 

5. Summary 

An epoch in low-mass star formation prior to 
the well-known Class stage is now well estab- 
lished. In this prior stage, which is now labeled 
a prestellar core, there is evidence for collapse to- 
wards a central condensation although this con- 
densation is still at the low temperature of the 
cloud. The observation of molecular column den- 
sities and radial distributions of these column den- 
sities in prestellar cores yields information on the 
nature of the collapse if it can be properly inter- 
preted. For such an interpretation, a detailed un- 
derstanding of the chemical processes is needed. 
In our first paper on the subject (Paper I), we 
modeled the gas-phase chemistry occurring dur- 
ing collapse according to the Larson-Penston so- 
lution. We found that both the average molec- 
ular column densities and their radial distribu- 
tions for CCS, CO, and N 2 H + in the prestellar 
core L1544 were only partially reproduced by the 
model. The most salient disagreement occurred 
in the column density of the molecular ion N 2 H + , 
which was determined to be much lower than ob- 
served. Slower collapse models did not improve 
the situation; N 2 H+ is enhanced but CCS column 
densities become significantly lower than observed. 
More recently, Li et al. (2002) performed similar 
calculations with two important additions: a small 
number of chemical reactions occurring on grain 
surfaces and a magnetic field to retard collapse. 
They found better agreement with observations 
in L1544, although it remained somewhat unclear 
whether the source of the improvement was the in- 
clusion of surface chemistry or the magnetic field, 
or both. 

In the present paper, we have included a com- 
plete set of surface chemical reactions into our 
chemical network and redone the calculations. 
There is now much more observational informa- 
tion with which to compare results - both in LI 544 
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and in other condensations. We find that our im- 
proved Larson-Penston model is now able to re- 
produce most of the observational data, especially 
in L1544. Furthermore, our current model results 
are in better agreement with Li et al. (2002) than 
shown by a comparison between Li et al. (2002) 
and our Paper I; ammonia and N2H + are more 
abundant, and the dependence of CCS column 
density on collapse time scale is smaller in the cur- 
rent model. The radial distributions of molecules 
are in reasonable agreement with those of Li et 
al. (2002) for their magnetized core, whether or 
not we use our standard case or slow the collapse 
down somewhat to mimic the magnetized result. 
Interestingly, the agreement is not as good with 
their unmagnctized core. Specifically, their ra- 
dial column density distributions of CS and CCS 
are strongly peaked at the core center in the non- 
magnetic case whereas ours are not. One must 
remember, in this comparison, that they use a sub- 
set of the UMIST network and we use the NSM 
network; differences in the abundances of sulfur- 
bearing species between the two networks are well 
known. On the whole, however, we may conclude 
that the inclusion of surface chemistry, however 
uncertain it may be, is clearly necessary 

How good is the current state of agreement be- 
tween our model and observational results? Re- 
stricting attention to L1544, we find that both 
the radial column density distributions of the well- 
studied species CS, CCS, N 2 H + and CO and aver- 
aged column densities of 9 species are in the main 
reasonably fit by a variety of models, including 
the standard Larson-Penston case and variations 
with a small degree of delayed collapse and/or re- 
duced sticking efficiencies of gas-phase molecules 
on grains. Specifically, our models can reproduce 
averaged column densities of up to eight species 
to within a factor of three. When the collapse 
is delayed by a factor of ten, however, the agree- 
ment is less good. The inferred radial distributions 
expressed in terms of fractional abundances are 
somewhat harder to reproduce with our standard 
initial conditions and an adsorption energy for CO 
based on the value for pure CO ice. A Moderately 
slow collapse model (/ ~ 3) in which at least some 
of the nitrogen is initially in its molecular form 
and the adsorption energy of CO is increased sig- 
nificantly to its value for pure water ice is much 
more successful. 



We have also attempted to consider a larger 
sample of prestellar cores and see how our mod- 
els can help to understand the varying degrees of 
evolution of the cores using two primary indica- 
tors of evolution: the central density of the core, 
and the central column density of the ion N 2 H + , 
which is found theoretically to increase monoton- 
ically with increasing core density. With these in- 
dicators, we are able to qualitatively reproduce the 
observed dependence of most of the molecular col- 
umn densities. The column density of NH 3 tends 
to increase, and those of CCS and CS to decrease 
with core evolution. The non-spherical geometry 
of cores should also be taken into account, how- 
ever, for more quantitative comparisons of species 
which are heavily depleted in the central regions. 

Another good indicator of evolution is pre- 
dicted to be the D/H ratio of molecular species, 
since fractionation is enhanced when density in- 
creases and heavy species are adsorbed onto grain 
surfaces. Indeed a plot of the observational 
data shows that the ratios of DNC/HNC and 
DCO+/HCO+ increase with the N 2 H+ column 
density. The dependence of DNC/HNC on the 
N 2 H + column density is in the main reproduced 
well by our standard Larson-Penston core and 
slight variations in collapse rate. The column 
density of N 2 D + and the high column density ra- 
tio of N 2 D+/N 2 H+ observed in L1544 are also 
reproduced in our models. There is one major 
problem though: the extended NSM network used 
in our calculations docs not contain some new 
rate coefficients measured by Gerlich, Herbst, & 
Roueff (2002); the effect of including these rates, 
which remain to be confirmed by other experi- 
ments, is to lower the fractionation to levels un- 
acceptable unless very slow collapse is assumed. 
Even so, the agreement is not as good as for the 
plain Larson-Penston case with the NSM network. 
This puzzling disagreement mimics that found by 
Roberts et al. (2002) for static cores. It is clear 
that deuterium fractionation is not yet completely 
understood. 

As more prestellar cores are observed in greater 
degrees of detail, it should be possible to distin- 
guish among chemical treatments better and to 
use model calculations to better constrain the evo- 
lutionary dynamics. Modeling of non-sphcrical 
cores and direct coupling of the large chemical 
network, (magneto-)hydrodynamic calculation of 
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cores and radiation transfer are desirable for com- 
parison with observed line profiles and intensity 
maps. Since surface chemistry appears to play 
an important role, improvements in the currently 
quite approximate treatment of surface processes 
would appear to be a necessity. Such improve- 
ments are currently being developed by several 
groups including our own. 
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Fig. 1. — Temporal variation of molecular abun- 
dances relative to hydrogen nuclei in a fluid cle- 
ment that migrates from 8.2 x 10 3 to 1.0 x 10 3 
AU in 2.00 x 10 5 yr. Only the later stages of the 
evolution are shown. The top panel (a) depicts 
the results of our current paper with grain-surface 
reactions, while the bottom panel (b) shows re- 
sults obtained without grain-surface reactions, as 
in Paper I. 



4000 S0O0 12000 4000 8000 12000 
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Fig. 2. — Distributions (n(i)/nn) of molecules at 
t = 1.52 x 10 5 yr, 1.89 x 10 5 yr, and 2.00 x 10 5 
yr (from top to bottom), when the central density 
of the core is 3 x 10 5 cm~ 3 , 3 x 10 6 cm~ 3 , and 
3 x 10 7 cm~ 3 , respectively. The left panels show 
distributions of gaseous molecules, while the right 
panels show those of adsorbed species ( "ices" ) on 
grain surfaces. 
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Fig. 3. — Molecular column densities as a func- 
tion of impact parameter (radial distance) from 
the core center when the central density of the 
core is 3 x 10 5 cm~ 3 (dotted lines), 3 x 10 6 cm~ 3 
(solid lines), and 3 x 10 7 cm~ 3 (dashed lines). 
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Fig. 4. — Molecular column densities as a func- 
tion of impact parameter when the central density 
of the core is 3 x 10 6 cm~ 3 . The solid lines rep- 
resent the Larson-Penston core, while the dashed 
and dotted lines represent cores collapsing more 
slowly by factors / of 3 and 10, respectively. The 
sticking probability is 1.0 for all thick lines, and 
1/f for thin lines. 
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Fig. 6. — Radial distributions of the abundances 
of deuterated and normal species in the Larson- 
Penston core when the central density of the core 
is 3 x 10 5 cm" 3 (top), 3 x 10 6 cm" 3 (middle), 
and 3 x 10 7 cm" 3 (bottom). The left panels show 
the abundances of gas-phase species, and the right 
panels show those of adsorbed species on grain sur- 
faces. 



Fig. 5. — Molecular column densities in the 
Larson-Penston core when the central density is 
3 x 10 6 cm" 3 . The initial central density of the 
core (cm -3 ) is 2 x 10 4 for the solid lines, 2 x 10 3 
for dashed lines, and 2 x 10 2 for dotted lines. The 
dashed-dotted lines show the molecular column 
densities obtained using the UMIST-based chemi- 
cal reaction network. 




Fig. 7. — Column density ratios of some deuter- 
ated isotopomers to normal species as a function 
of impact parameter from the core center when the 
central density is 3 x 10 6 cm -3 . Solid lines rep- 
resent the Larson-Penston core, while dashed and 
dotted lines represent cores collapsing more slowly 
by factors / of 3 and 10, respectively. The stick- 
ing probability is 1.0 for thick lines, and 1/f for 
thin lines. Dashed-dotted lines show the molecular 
column densities in the Larson-Penston core, but 
with the "New Rates" for deuterium fractionation 
reactions. 
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Fig. 8. — Radial distribution of molecular abun- 
dances in L1544 estimated by Tafalla et al. (2002) 
(a), and obtained in our models when the central 
density of the core is nn = 3 x 10 6 cm -3 (b-f). 
Part (b) shows results from our Larson-Penston 
core, while the slow collapse model (/ = 3 and 
S = 1) is shown in (c). Parts (d-f) show results 
from revised models with nitrogen initially in its 
molecular form and a higher adsorption energy for 
CO. 
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Fig. 9. — Computed molecular column densities 
of assorted species as a function of the central 
density of the core. Column densities are surface- 
averaged over a radius of 4340 AU (which corre- 
sponds to 31" assuming the distance to the core 
to be 140 pc), 2000 AU, 3780 AU (27"), 5600 AU 
(40"), 2380 AU (17"), and 2000 AU from the core 
center for CCS, CO, N 2 H+, NH 3 , CS, and HCO+, 
respectively. The solid lines represent the Larson- 
Penston core, while the dashed and dotted lines 
represent more slowly collapsing cores by a factor 
/ of 3 and 10, respectively. The sticking probabil- 
ity is 1.0 for thick lines, and l/f for thin lines. 
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Fig. 10. — Column density ratios of deuterated to 
normal species as a function of suitably averaged 
N2H + column density. The solid lines represent 
the Larson-Penston core, while the dashed and 
dotted lines represent more slowly collapsing cores 
by a factor / of 3 and 10, respectively. The stick- 
ing probability is 1.0 for thick lines, and 1// for 
thin lines. The observational results are depicted 
with error bars; horizontal and vertical lines in- 
dicate uncertainties in the N2H + column density 
and molecular D/H ratio, respectively. 
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